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Abstract 



This document contains refereed supplementary information for the paper "Rapid planetes- 
imal formation in turbulent circumstellar discs". It contains 15 sections (§1.1 - §1-15) that 
address a number of subjects related to the main paper. Some of the subjects are highlighted 
here in the abstract. We describe in detail the Poisson solver used to find the self-potential 
of the solid particles, including a linear and a non-linear test problem (§1.3). Dissipative col- 
lisions remove energy from the motion of the particles by collisional cooling (§1.4), an effect 
that allows gravitational collapse to occur in somewhat less massive discs (§1.7). A resolution 
study of the gravitational collapse of the boulders is presented in §1.6. We find that gravita- 
tional collapse can occur in progressively less massive discs as the grid resolution is increased, 
likely due to the decreased smoothing of the particle-mesh self-gravity solver with increasing 
resolution. In §1.10 we show that it is in good agreement with the Goldreich & Ward (1973) 
stability analysis to form several-hundred-km-sized bodies, when the analysis is applied to 5 
AU and to regions of increased boulder column density. §11 is devoted to the measurement 
of random speeds and collision speeds between boulders. We find good agreement between 
our measurements and analytical theory for the random speeds, but the measured collision 
speeds are 3 times lower than expected from analytical theory. Higher resolution studies, and 
an improved analytical theory of collision speeds that takes into account epicyclic motion, 
will be needed to determine whether collision speeds have converged. In §1.12 we present 
models with no magnetic fields. The boulder layer still exhibits strong clumping, due to the 
streaming instability, if the global solids-to-gas ratio is increased by a factor 3. Gravitational 
collapse occurs as readily as in magnetised discs. 

Authors: 

Anders Johansen 1 , Jeffrey S. Oishi 2 ' 3 , Mordecai-Mark Mac Low 2,1 , Hubert Klahr 1 , Thomas 
Henning 1 & Andrew Youdin 4 

Affiliations: 

1. Max-Plan ck-Institut fur Astronomie, Konigstuhl 17, D-69117 Heidelberg, Germany 

2. Department of Astrophysics, American Museum of Natural History, 79th Street at Cen- 
tral Park West, New York, NY 10024-5192, USA 

3. also Department of Astronomy, University of Virginia, Charlottesville, VA, USA 

4. Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George 
Street, Toronto, Ontario M5S 3H8, Canada 



1 Supplementary Discussion 



1.1 Pencil Code 

The Pencil CodeP is a finite difference code that uses symmetric derivatives of sixth order 
in space and a third order Runge-Kutta time-stepping scheme. The difference equations as 
implemented are formally dissipation free, having phase errors but no amplitude error. Only 
a small amount of numerical dissipation is introduced from time-stepping the advection term. 
Therefore, one must explicitly add dissipation to the dynamical equations to suppress nu- 
merically unstable modes near the grid scale and to dissipate the turbulent energy that is 
(in our case) released from the Keplerian shear by the Reynolds and Maxwell stresses. For 
this purpose, we use sixth order hyperdiffusivity operators, where the usual V 2 diffusivity 
operator is replaced with a V 6 operator! 32 * 33 !. Hyperdiffusivity dissipates energy at high wave 
numbers - at the smallest scales in the simulation - but preserves energy at low wave numbers. 
Hyperviscosity and magnetic hyperresistivity have been used extensively to study the proper- 
ties of forced magnetohydrodynamic turbulence (see Brandenburg & Sarson^ and references 
therein). They are designed to affect large scales as little as possible by dissipation, thus 
widening the inertial range beyond what can be achieved with a regular viscosity operator 
while still maintaining numerical stability. We have tested that particle overdensities occur 
for both the stringent hyperviscosity scheme used by Haugen & Brandenburg^ and for the 
simplified scheme described in Johansen & Klalrr^. 

Possible side effects of using hyperviscosity and hyperresistivity include an artificial increase 
in the bottleneck effect^, a physical effect in turbulence where energy piles up around the 
dissipative scale, and a higher saturation level for dynamo-generated magnetic fields in helical 
flows compared to what is seen when using a regular viscosity operator^. The bottleneck 
effect is unlikely to be relevant to the dynamics of boulders in turbulence, since marginally 
coupled particles are mostly affected by turbulent structures at the largest scales of our 
simulation where the bottleneck effect is unimportant. The saturation level of turbulence 
driven by the magnetorotational instability (MRI) is indeed affected by the numerical scheme 
and by the dissipation scheme, but the Pencil Code agrees well with other grid codes regarding 
the statistical properties of MRI turbulenc d 14 * 36 l 

1.2 Drag force 

1.2.1 Drag force algorithm 

The computation of drag forces between Lagrangian particles and an Eulerian grid requires 
some care to avoid spurious accelerations and to ensure momentum conservation. Small 
errors in the gas velocity can be dangerously amplified by the subtraction of highly correlated 
particle velocities. Tests of our drag force algorithm are described in detail elsewhere^. It 
involves three steps: 

1. Interpolate gas velocities at particle positions. 
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2. Calculate the drag force on the particles from the gas in nearby cells. 



3. Assign the back-reaction force to the gas from particles in nearby cells. 

For the first step, interpolation, we begin with gas velocities, defined on a uniform grid 
where the index j labels the cells centred on positions We interpolate to the particle 

positions, x®, using a weight function, Wi, as 



u (x«) = Wi(x® - cc (5 V W) • (1) 

3 

The weight function is normalised as J2j Wi(x^ — x^') = 1, for any x®, and has non-zero 
contributions only from the cells in the immediate vicinity of 

The second step, calculating the drag acceleration on particle i, 
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/« = t,W-«(x«) , (2) 

1 Tf L J 

is trivial once the relevant quantities are defined. We assume that the drag force is propor- 
tional to the velocity difference between gas and particles, i.e. that the friction time Tf is 
independent of the velocity difference. 

Finally, we calculate the back-reaction drag force, fg, on the gas in cell j. Assigning the 
particle velocities to a grid risks violating momentum conservation. Instead we use Newton's 
third law to directly assign the force on the particles back to the gas, 

= E w ^ - . (3) 

Pg V ce u i 

(i) 

where m p is the mass of a superparticle, p g is the gas density in cell j, and V^ e ii is the 
volume of a grid cell. The assignment function Wa obeys the same conditions as Wj, so that 
only particles in a given cell or in that cell's nearby neighbours contribute to the sum. We 
opt for the second order Triangular Shaped Cloud (TSC) assignment scheme®, which uses 
an identical function for assignment and for interpolation, Wa = W\. The TSC spreads the 
influence of particles and grid points to three grid points in each direction, for a total of 27 
points in three-dimensional (3-D) simulations. 

The interpolation errors associated with the TSC assignment is found by considering a periodic 
function (of arbitrary phase) sampled at the grid points®^. The result is that the assigned 
amplitude of a single Fourier component at scale k relative to the actual amplitude is to 
second order 1 — (Ak) 2 /8, where A is the linear size of a grid cell. Thus already at 5 grid 
cells there is significant smoothing. This smoothing is found to have an influence on the 



gravitational collapse, especially at crude resolution, see £1.6 



The time-step constraint set by the drag force is <5idrag = Tf/(1 + e), where e is the local 
solids-to-gas ratio^. As the solids-to-gas ratio increases, the allowed time-step decreases. We 
have made sure that the friction time-step never dominates over the Courant time-step of 
the code by artificially increasing the friction time in regions of high solids-to-gas ratio (e > 
100). We have experimented with the threshold and found no improved collapse for a higher 
threshold, presumably because inelastic collisions dominates the kinetic energy dissipation at 
these densities anyway. 
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1.2.2 Particle sizes and drag regimes 



In our work, as in most theoretical work, we characterise particles in terms of their dimen- 
sionless friction time i?K r f (where i?K is the Keplerian orbital frequency). The translation 
to a particle size depends on the assumed disc model. The smallest particles are subject to 
Epstein drag, valid when the particle radius is smaller than the mean free path of the gas 
(see below), with 

^ kt (Ep) = %» = ^P^a (4) 

where the second step applies in the disc midplane. Here p 9 is the material density of the 
solids, a is the radius of a solid body, p g is the gas density, c s is the sound speed, while 
E g is the column density of gas. Epstein drag depends on gas density, but in practice the 
gas density fluctuations are negligible (order 1%) in our subsonic flows, so we ignore them. 
Solving for particle size gives 

V^p. ~ 3 ° Cm ^ KT f U50gcm-2j Ugcm-3j [ 5 AU J > W 

where the normalisation of the gas surface density at 5 AU, S g 5, and the power law slope 
follows the minimum mass solar nebula modeF^. Applying a simulation with given Tf values 
to different disc radii changes the relevant particle sizes. Considering, as in the main text, 
r = 5AU and S g = 300gcm~ 2 yields particle sizes of a = 60,45,30,15 cm, for i?K r f = 
1.0,0.75,0.5,0.25, respectively. When the model is applied to the outer solar nebula at r = 
40 AU instead, the corresponding particle sizes are as low as a few centimetres (per J?K r f)- 

The Epstein regime of free molecular (or Knudsen) flow ceases to apply once the particle 
radius exceeds (9/4 of) the gas mean free path, 

A = _Jf_ = VpM^ (6) 

PgCmol ^g^mol 

v^Ogcm-V 0.04 UAUy ' { ' 

where p. = 3.9 x 10 -24 g is the mean molecular weight and <r mo i = 2 x 10 _15 cm 2 is the molecular 
cross section of molecular hydrogen ^ 39 ! 40 ! . The radial dependence does not include flaring of 
the aspect ratio, H/r. Epstein drag applies as long as 

(Ep) 9. / p. \ / E gt5 X- 2 / H/r\ (jr^y 

K f 2 r 2 (j mol Ugcm-sJ ^I50gcm- 2 y ^0.04/ \5 AU J ' 1 ' 

Thus Epstein drag is the relevant regime for the application of our model to 5 AU. 

Extrapolation of our model to higher column densities (or regions closer to the star) requires 
consideration of Stokes drag, for which 

n K7 ft) = fl (E P )4a = 4 P .aV mol 
K f K f 9A 9pH V ; 
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Since Stokes drag is also linear in the relative velocity between gas and solids, our model 
applies with a different scaling of particle size, 



9J? K r f (St V# 
4/}.<r mo i 



,1/2 



80 cm 



(St) 



P* 



2 g cm 



0.04/ V5AU 



.1/2 



(10) 



Actually since Stokes drag is independent of gas density (including fluctuations), our model 
is more exact in this regime, but modelling particles in the Stokes regime would require a 



different treatment of collisional cooling (see £1.4). 



1.2.3 Energy dissipation by drag force 

The dissipation of kinetic energy caused by drag force can be calculated as 

V ot J drag V ot J drag 

where u and w are the respective gas and particle velocity fields. We refer to references^ 
anf jH2] f or considerations of the effect of drag force damping on the dynamics of particles. 
Inserting the fluid expressions for the drag force acceleration, 

^ = -^(u-w), (12) 
at Tf 

dw 1 . . , 

8T = -r, [m - U) < < 13 ' 

yields 

e kin = -^\w-u\ 2 . (14) 

Dissipation of kinetic energy by drag forces comes automatically when applying the momentum- 
conserving drag force. We usually assume that the dissipated energy can immediately radiate 
away efficiently, keeping the temperature of the gas constant. We show, however, in §1.13 that 



even if all the released energy remains locally as heat, the resulting temperature increase of 
the gas is insignificant and has no influence on the initial stages of the gravitational collapse. 

As long as there is a velocity difference between the solid particles and the surrounding gas, 



then drag force cooling is efficient according to equation (14). But the relative speed of gas 
and particles in drag force equilibrium decreases with increasing solids-to-gas ratio^, as the 
solid particles entrain the gas. In that case inelastic collisions take over as the dominant 



cooling mechanism (see £1.4) 



1.3 Self-gravity solver 



We compute the gravitational potential of the particles by determining a particle density on 
the mesh and solving the Poisson equation for this assigned density field. The gravitational 
acceleration is then interpolated back to the positions of the particles. The particles are 
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assigned to the mesh using the TSC scheme^, spreading each particle over the 27 nearest 
grid points (in 3-D), and the gravitational acceleration is added back to the particles using a 
second order spline interpolation to avoid any risk of self- acceleration of the particles. 

The Poisson equation, 

V 2 <2> = 4vrG/j p , (15) 

where G is the gravitational constant and p p is the assigned particle density field, is solved 
using a Fourier method. The solution for a single Fourier density mode of wave number 
k = (k x ,k y ,k z ) and complex Fourier amplitude pu is 

$k = p— (16) 

for k = \k\ > 0. The full potential is then <P(x) = X)fc^fcexp[ifc • x\. 

The radial x-direction in the shearing sheet is not strictly periodic, but is rather shear periodic 
with the connected points at the inner and outer boundary separated by the distance Ay(t) = 
mod[(3/2)/?K-^a;i) Ly] in the y-direction. We follow the algorithm of Gammie 20 to include 
shear-periodic boundaries in the Fourier method, with the difference that we perform the 
necessary interpolation in Fourier space rather than in real space, as suggested by C. McNally 
at |http : / /i mp . mcmaste rTca/~colinm/ ism/rotf f t . html| First we apply a discrete Fast 
Fourier Transform (FFT) in the periodic y-direction. We then shift the entire y-direction 
by the amount 5y(x) = Ay(t)x/L x to make the x-direction periodic, before proceeding with 
discrete FFTs along x and then z. Performing the 5y(x) shift in Fourier space (essentially 
using a Fourier interpolation method) has the advantage over polynomial interpolation that it 
is continuous and smooth in all its derivatives. After solving the Poisson equation in Fourier 



space, using the algebraic solution from equation (16), we transform the potential back to 
real space in the opposite order, shifting the y-direction back to the sheared frame on the 
way. 



The gravitational constant G in equation (15) must be defined in code units. We adopt the 



unit system i?K = c s = H = p gjm id = 1, where the parameters are respectively the Keplerian 
frequency, isothermal sound speed, pressure scale height, and midplane gas density of the 
disc. The non-dimensional form of the Poisson equation is 

{HVfZ/cl = G—^— , (17) 

Pg,mid 

where 

G = ATiGp gMA Q K 2 (18) 

is a self-gravity parameter which relates to the Toomre parameter^ Q of self-gravitating gas 
discs as 

Qml&CT 1 . (19) 

In case of vertical hydrostatic equilibrium the connection between the midplane density p^rmd 
and the column density U g is p g , m id = S S /(V^H). Introducing further the dimensionless 
scale-height-to-radius parameter H/r yields the power law connection 
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Here we have made use of the identity fi\ = GM^/r 3 , where M* is the mass of the central 
gravity source and r is the orbital radius to a considered point in the disc. 



1.3.1 Testing the Poisson solver 

We validate here the FFT Poisson solver for a complex test case of a self-gravitating, shearing 
particle wave with gas drag. This is not a standard test problem, but we found it very useful 
for validating the particle-mesh self-gravity scheme of the Pencil Code. The particle density, 
velocity and self-potential are written as p p = po + p' p , w = wq + w' , <P = <&q + Here 

wq = Uy y = —(3/2)f2xxy is the Keplerian shear flow. The linearised continuity equation, 
equation of motion and Poisson equation are 

d t p' p + u^d y p' p = -p (d x w' x + d y w' y ), (21) 

d t w' x + ufd y w' x = 2f2 K w' y - 8 X $' - -w' x , (22) 

dtW'y+U^dyW'y = -\{2 K w' x -dy& (23) 

Z Tf 

V 2 <2>' = 4vrGpp. (24) 

Here we have assumed particle motion to take place in the radial-azimuthal plane and zero 
gas velocity, u x = u y = 0, in the friction terms. Because of the Keplerian shear flow, the 
coefficients of the linearised equations are not independent oft, and thus an eigenmode analysis 
is not possible. Instead, we assume a separable solution q(t,x,y) = q(t) exp[i(k x (t)x + k y y)\ 
for each dynamical variable. The time derivative of q(t, x, y) is 

q(t, x, y) = [d t q(t) + q(t)ixd t k x (t)} exp[i(k x (t)x + k y y)\ . (25) 



By setting xdtk x (t) = —Uy°^k y , we cancel the u y °^d y q terms in equations (21)-(23), leaving a 
system of equations for (p p ,w x ,w y ), 



dt 



-p i(k x w x + kyWy) , (26) 
dwx nn + . 47riGk x pp 1 „ 

-df = ™ KW y + i^Tkj ~ n Wx ' (27) 

dwy 1 A AmGkyPp 1 



;Vkw x + , 2 , ^ p - -w y , (28) 



together with k x (t) = k x (to) + (3/2) filthy. 

We solve this system of ordinary differential equations numerically using a third order Runge- 
Kutta method to follow the temporal evolution of a non-axisymmetric wave with the initial 
condition k x = —1, k y = l,w x = w y = 0, p p = 10~ 4 , and G = i?K = Po = r f = 1- This semi- 
analytic solution is then compared to the evolution obtained with the full solver of the Pencil 
Code (64 2 grid points and 4 particles per cell). Because we use particles to represent the solids 
and a grid-based FFT solver to calculate the gravitational potential, the comparison provides 
an excellent of test of all physics and numerical schemes relevant to our investigation: nu- 
merical particles, grid-based FFT self-gravity, shear and shear-periodic boundary conditions. 
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Figure 1: The evolution of a self-gravitating shear wave of solid particles. The plot shows 
a comparison between the semi-analytic solution to the linearised equation system (grey 
lines) and the solution obtained with the full solver of the Pencil Code (dark dotted, dashed 
and dot-dashed lines) for the amplitude of the particle density p p and of the particle velocity 
components w x and w y . There is excellent agreement between the numerical and the analytical 
solutions in the entire linear range. After p p ~ 0.1 the analytical solution is no longer valid 
and the two solutions diverge. 



We show the evolution of the analytical solution and of the solution obtained with the Pencil 
Code in Fig. [TJ There is an excellent agreement between the two solutions up to around 
Pp = 0.1 where the problem enters the nonlinear regime and the analytical solution loses its 
validity (but the Pencil Code solution does not). 

Next we derive an analytical solution to a fully non-linear one-dimensional gravitational col- 
lapse problem, following Spitzer^, to test the Poisson solver outside the linear regime. 

Consider the equation system governing a self-gravitating gas in an infinitely extended one- 
dimensional space along the z-axis: 

du z du z d<P ndlnp . , 

<91np dlnp du z 

~dT + u ^ = ~W m 

d 2 <P 

= 4irGp. (31) 



dz 2 

We search for stationary solutions with u z = p = u z = 0. This yields a second order 
differential equation in p, 

9 d 2 In p „ , . 

c 2 s ^-f+4irGp = 0. 32 
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Figure 2: Comparison of numerical and analytical solution to the Spitzer^ one-dimensional, 
non-linear collapse problem. The agreement is very good in high density regions, but there is 
disagreement in the underdense regions because the analytical solution assumes infinite space 
rather than a periodic domain. 

The general solution is 



is given by the column density of the gas and zq is an arbitrary constant. This is an exact 
solution to the full non-linear equation system and thus complements the linear test problem 
described in the previous section. We let the Pencil Code start with a Jeans unstable mode 
of wavenumber k z = 0.5 in a periodic z-space of length L z = 4ir covered by 64 grid points (we 
used shock viscosity to keep the code stable during the initial collapse). We show in Fig. [2] 
a comparison between the equilibrium state found by direct time integration of the equation 
system by the Pencil Code and the analytical equilibrium solution. There is an excellent 
agreement between the two solutions in regions of high and moderate density, whereas the 
underdense regions contain too much mass in the Pencil Code solution. This is however just 
an artefact of the comparison: the analytical solution works in infinitely extended space, 
whereas we assumed periodic boundary conditions in the numerical solution. 




(33) 



where 




(34) 
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1.4 Collisional cooling 



In the models that include collisional cooling we let unresolved collisions damp the velocity 
dispersion of the particles within each grid cell. The friction time in the Epstein^ regime is 

a »P» ar\ 

r f = , (35) 

where we follow the nomenclature of equation Q. The collisional time-scale is assumed to 
follow a simple scaling with the Epstein friction time, 

PgC s Tf = PpCpT co ll , (36) 

where p p is the bulk mass density of solids and c p is the velocity dispersion of the boulders in a 
grid cell. We ignore for simplicity any proportionality factors of order unity (for the complete 
analytical expression of the cooling time-scale, see Garzo & DuftyBBJ). If particles are in the 



Stokes drag regime rather than Epstein (see { 1.2.2), then the collisional cooling time-scale 



actually decreases relative to the friction time, but for simplicity we do not model this regime 



here. Equation (36) allows us to calculate the collisional time-scale from the friction time as 



^_ *f (37) 

C ° U (Cp/c s )( / 9 p /p g ) ' 

without any reference to the radius and solid density of the particles. We discuss below how 
this expression changes when particles of several sizes collide. 

Taking a velocity dispersion of c p /c s ~ 1 x 10~ 2 for the random motion of particles within 
a grid cell gives a threshold solids-to-gas ratio of e ~ 100 where collisions between boulders 
becomes as important as collisions between boulders and gas molecules. The collisional cooling 
is implemented as a simple term that reduces the velocity of a particle relative to the mean 
velocity of the particles in a grid cell on a collisional time-scale, 

9 ^ = - 1 ^(v^-v^), (38) 

Ot T C oll 

where is the average particle velocity in cell k in which particle i is situated. The 
coefficient of restitution C res is a measure of the degree of energy conservation in the impacts. 
We take C res = 0.1, corresponding to a loss of 90% of the relative velocity in each collision, 
a reasonable value for collisions between macroscopic solids^, giving a collisional cooling 
time-scale that is approximately the same as the collisional time-scale 



For collisions between particles of multiple sizes we replace equation (37) with the equation 



Tcol j = f (39) 

(cp/cs) E;W/>g) ' 

where pi is the local bulk density of solids of species i. The size-modified friction time r f * is 
an average of the friction times of the individual species weighted with the local bulk density 
of the individual species, 

r * = . (40) 
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Using Tf as the effective friction time, we can define a single cooling time-scale for all particle 
size bins. This ensures that the parameterised collisional cooling conserves momentum. For 
a single particle size the effective friction time reduces to the friction time, r f * = Tf. 

In the numerical simulations we calculate the collision speed of particles in a given grid cell i 
by first calculating the mean velocity in the cell tfW. We then calculate the average particle 
collision speed as 

cf = |«« -«W| . (41) 



This value is then plugged into equation (37) to get the collisional time-scale. This approach 
must be taken when using superparticles, as collisions cannot directly be followed as could, 
in principle, be done for individual particles. 



1.5 Initial condition 



Our simulations are initialised with a uniform density gas, p g = 1, with isothermal sound 
speed c s = 1. The superparticles are randomly distributed in space, and the mass of each 
superparticle is set such that the column density of solids is £ p = O.OlV^-Hpg.mid- Here we 
have assumed that prior to the start of the simulation, sufficient sedimentation has occurred 
to bring all the solids within the computational domain. The particles are given random 
velocities with a Gaussian distribution of width Sv = 1 x 10 _3 c s . The domain is cubic with 
L x = L y = L z = 1.32H, where H is the pressure scale height of the disc. We use periodic 
boundary conditions in the azimuthal y-direction and in the vertical z-direction and shearing 
periodic boundary conditions in the radial x-directionS^l. The Keplerian rotation frequency 
at the central radius of the grid is J?k = 1- Assuming an orbital distance of r = 5AU, a 
gas column density of U s = 300gcm~ 2 and a scale-height-to-radius ratio of H/r = 0.04 then 
gives H = 3 x 10 12 cm, p g = 4 x 10~ n gcm -3 and c s = 500 ms" 1 as physical properties of 
the disc. The column density of solids is U p = 3gcm~ 2 for a global solids-to-gas ratio of 
e = 0.01. 



1.5.1 Radial pressure support 



The assumed radial pressure gradient, din P/ din r, causes radial drift of the boulders follow- 
ing the expression^ 



(d In P/d In r){H/r) 



2Av 



f2 K T { + {Q K Ti)- 1 ~ H ^KTf + (^KTf)- 1 

in the absence of collective drag force effects on the gas (i.e. vanishing solids-to-gas ratio). 
With {d In P/d In r){H/r) = —0.04 this gives v x = Av = — 0.02c s for marginally coupled 
boulders with J?K T f = 1- We show in fjl.14 the effect of less and more radial pressure support. 



1.5.2 Magnetic fields 

The Pencil Code solves the induction equation for the magnetic vector potential A, keeping 
the magnetic field B = V x A solenoidal (i.e. V • B = 0) per construction. We initialise 
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magnetic fields by setting the initial vector potential 



/2tt \ 2tt \ /2vr , 

A cos ( —x J cos I — y j cos ^— z ) y, ( K!) 



with Ao = 0.04, which leads to an initial zero-net flux field with a maximum thermal to 
magnetic pressure ratio f3 ~ 55. For the 256 3 run we add a very weak external vertical 
magnetic field of constant /3 ext = 2 x 10 5 . We found this latter field to be necessary to 
maintain a constant turbulent viscosity of a = 10~ 3 throughout the simulation. 



1.5.3 Disc mass 



We set for the 256 3 self- gravity simulation presented in the main text a value of G = 0.1, 
corresponding to a disc with Toomre^l parameter Q ~ 16, an order of magnitude below 
the limit for global gravitational instability of the gas at G = 1 (at very large scales in an 
infinitely extended disc). For this reason we do not let the gas feel self-gravity, nor do we 
let the gas contribute to the gravitational potential in which the particles move. We have 
performed simulations that included the self-gravity of the gas, but found no significant gas 
overdensity in the gravitationally contracting particle clusters, due to the strong pressure 
support of the gas. Our choice of G gives a gas surface density of S g = 300 gcm~ 2 at 
r = 5AU (see eq. 20), approximately two times the minimum mass solar nebula (MMSJNpSl). 
With moderate radial pressure support Av = — 0.05c s the column density limit is twice as 



high as for Av = — 0.02c s , see [1.14 These are not unrealistic values for the actual column 
density of gas in the solar nebula, since the MMSN is a very conservative estimate of the 
column density that is calculated by adding up the heavy elements in the present Solar 
System and then dividing by the assumed interstellar solids-to-gas ratio. Any loss of material 
that occurred due to radial drift will then lead to an underestimate of the actual column 
density. Radial drift can also enhance the solids-to-gas ratio in the inner 10 AU by as much 
as an order of magnitude™^. 



1.5.4 Particle sizes 



We assume in all the models described here a global solids-to-gas ratio of eq = 0.01. There 
may be twice as much solid material present due to the condensation of volatiles into ice 
beyond the snow line at a few AU from the protostar. Coagulation models generally yield 
a particle mass distribution that spans two orders of magnitude in radius'^SEH. Our particle 
size distribution has size bin boundaries spanning an order of magnitude from l?K r f = 0.125 
to f^KTf = 1.125. The assumption is then that the simulated boulders represent the upper 
half of the actual size range and that the lower half is present in smaller particles that are 
not effectively concentrated by transient high pressures and streaming instability and also are 
unable to participate in the collapse due to their strong coupling to the gas. 
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Figure 3: 1-D spectra of the gas velocity components with three different resolutions along 
the rows, the three different velocity components along the columns, and spectra taken in 
three different directions shown in each panel, in simulations with no particles. The different 
colours show power along the three coordinate directions, while the dotted line indicates a 
/j- 1 / 3 Kolmogorov law. 



1.6 Resolution study 
1.6.1 Power spectra 

Magnetorotational turbulence is known to produce a Kolmogorov-like power spectrum® with 
velocity amplitude on scale k going as oc A; -1 / 3 in the inertial range where numerical 
viscosity is insignificant. Increasing resolution leads to an increase in the inertial range. 
We show in Fig. [3] the modulus of the velocity amplitude as a function of wavenumber k 
normalised with the wavenumber of the largest scale of the box. A Kolmogorov A; -1 / 3 power 
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Table 1: Resolution study for fixed G. Col. (1): Mesh resolution. Col. (2): Number of 
superparticles in millions. Col. (3): Number of orbits with self-gravity. Col. (4): Measured 
turbulent viscosity. Col. (6): Self-gravity parameter. Col. (7): Corresponding Toomre Q ~ 
1.6G -1 . Col. (8): Number of clusters at the end of the simulation. Col. (9): Accretion rate 
of the most massive cluster in Ceres masses per orbit. Col. (10): Accretion rate normalised 
with Go.i = G/0.1. Col. (11): Total mass in the four most massive gravitationally bound 
clusters at 7 orbits after self-gravity is turned on. 



law is indicated with a dotted line for reference. The turbulence is notably weaker in the 
z-direction than in the x- and y-directions. It also clear from Fig. [3] that the inertial range 
expands as expected with increasing resolution. Resolving the turbulent scales is important 
for resolving the local velocity dispersion ("temperature") of the particles and thus for the 
gravitational collapse (see £1.11). 



1.6.2 Accretion 



We check the convergence of the accreted mass using runs described in Table [TJ We compare 
runs at 64 3 and 128 3 zone resolution, using a column density (parameterised by the self- 
gravity parameter G) that allows the 64 3 simulation to undergo gravitational collapse into 
discrete clusters. We find that the 128 3 run forms four clusters while the 64 3 run only forms 
two. However, the most massive cluster accretes faster in the low resolution case, because 
it does not have to compete for particles with other clusters. In Fig. [4] we plot the mass in 
the four most massive clusters as a function of time. The clusters in the 128 3 run accrete 
at approximately 50% higher rate than in the 64 3 run, because there are more clusters upon 
which particles can accrete. This is also evident from the last column of Table [T] where we 
write the total mass of the four most massive, bound clusters at a time of 7 orbits after self- 
gravity is turned on. In Fig. [5] we plot the maximum bulk density of particles and the mass 
of the most massive bound cluster as a function of time (similar to Fig. 3 of the main text). 
It is clear that the most massive cluster contains more mass and accretes faster in the 64 3 
simulation. At the same time the maximum density is higher in the 128 3 simulation, because 



the increased resolution leads to less smoothing in the particle-mesh scheme (see {1.2). 



We emphasise that the gravitational collapse of overdense seeds in turbulence is inherently a 
stochastic process, because of the temporal variability of the overdensities, so that convergence 
is not expected in the strict meaning of the term. Only convergence to mean values can be 
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Figure 4: Total mass in the four most massive gravitationally bound clusters, in units of 
Ceres masses, plotted versus time since self-gravity was turned on, for two resolutions but 
a fixed self- gravity parameter G that allows gravitational collapse at both 64 3 and at 128 3 . 
The higher resolution run accretes at a 50% higher rate onto four clusters condensed out of 
the midplane layer, while only two clusters formed in the 64 3 run. 

expected over a large ensemble of simulations. The limited number of tests that we have 
been able to perform within available computational resources should therefore be taken as 
qualitatively suggestive rather than conclusive. 

1.6.3 Gravitational collapse 

We have also done a resolution study of the column density required for gravitational col- 
lapse of the particle layer to occur, with results given in Table [2j The midplane layer in all 
runs presented here is resolved with an initial average of around four particles per grid cell, 
enough for the TSC scheme to define meaningful particle densities for the self-gravity solver. 
Collapsing regions have many more particles in each cell, of course. The column density limit 
for forming gravitationally bound clusters (see column 6 of Table [2]) decreases with increasing 
grid resolution at all resolutions we have been able to study to date, as determined by varying 
the column density in increments of AZ" g = 300 g cm~ 3 , equivalent to AG = 0.1. 

The observed decrease in limiting column density appears to be partly because of a decrease in 
the strength of the magnetorotational turbulence with increasing resolution, at least at lower 
resolution. The strength of the turbulence can be parameterised by the effective turbulent 
viscosity, which drops from a = 0.002 at 64 3 to a = 0.001 at 128 3 and 256 3 , allowing a slightly 
denser midplane layer to form. 

A more important factor, though, is that the particle density field becomes better resolved on 
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Figure 5: The maximum bulk solid density and the mass of the most massive bound cluster 
as a function of time for a fixed value G = 0.4 at 64 3 (left, 1.25 x 10 5 particles) and at 128 3 
(right, 1 x 10 6 particles). The 64 3 run produces only 2 clusters (see Table [l]), but the most 
massive one accretes lots of solid material (7.5 Mc er es per orbit), whereas at the increased 
resolution of 128 3 four clusters form and they must compete for particles to accrete, thus the 
accretion rate of the most massive cluster decreases. The maximum bulk density of particles 
increases for higher resolution because the particle density field is better resolved by the TSC 
assignment scheme (see {1.2). Fig. [4] shows that the total mass accretion onto bound clusters 
in the 128 3 simulation is actually twice as high as in the 64 3 simulation. 



the grid at higher resolution. The TSC assignment scheme underestimates the amplitude of 
density modes with wavelength near the grid size^, so the self-gravity solver underestimates 
the gravitational acceleration caused by structures at those small scales. This can be seen in 
Fig. [6j where we show the maximum particle density produced by streaming instability and 
concentration in transient high pressures in models with neither self-gravity nor collisional 
cooling. Increasing linear resolution by a factor of two increases maximum density by around 
a factor four. Thus only higher resolution models reach the Toomre criterion for collapse at 
low column densities. 

We write in columns 12 and 13 of Table [2] the mass in the four most massive, bound clusters 
and the same parameter divided by Go.i = G/0.1. There is less (normalised) mass accretion 
at 256 3 than at lower resolutions. This is readily explained because collapse can take place 
in smaller regions when the resolution is increased. The clusters thus contain less mass and 
accrete slower, at least to begin with, although the 256 3 has not run long enough to reveal 
the behaviour over longer times. 

Another important effect is the resolution of the initial gravitational contraction. Fig. 2 in 
the main text shows that self-gravity results first in a radial contraction of the overdense seed 
bands. This contraction continues until the bands reach the Hill density where self-gravity 
dominates over tidal force, allowing a full non-axisymmetric collapse into bound clusters oc- 
curs. However, the self-gravity solver cannot follow collapse below a few grid cells, setting a 
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Table 2: Resolution study for gravitational collapse. Col. (1): Mesh resolution. Col. (2): 
Number of superparticles in millions. Col. (3): Number of orbits with self-gravity. Col. (4): 
Measured turbulent viscosity. Col. (6): Minimum self-gravity parameter where gravitationally 
bound clusters form (MMSN has G ~ 0.05 at r = 5 AU). Col. (7): Corresponding Toomre 
Q ~ 1.6G . Col. (8): Minimum self-gravity parameter where gravitationally bound clusters 
form without collisional cooling. Col. (9): Number of clusters at the end of the simulation. 
Col. (10): Accretion rate of the most massive cluster in Ceres masses per orbit. Col. (11): 
Accretion rate normalised with Gq.i = G/0.1. Col. (12): Total mass in the four most massive 
gravitationally bound clusters at 7 orbits after self-gravity is turned on. Col. (13): Total 
cluster mass normalised with Gq.i = G/0.1. 



resolution-dependent minimum width on the radial bands. Increasing the resolution decreases 
the final width of the bands, increasing their density, and thus allowing full collapse to occur 
in models with lower initial column densities. After collapse sets in, higher resolution al- 
lows collapse to higher densities, again resulting in higher peak densities in higher resolution 
models. 

Poisson noise due to the discrete nature of the superparticles on the other hand, appears to 
be quite insignificant for gravitational collapse. We have done tests with particle numbers as 
high as thirty particles per cell (in the mid-plane layer) in the 64 3 runs (10 6 particles) and 
eight particles per cell in the 128 3 runs (2 x 10 6 particles) and found the column density for 
collapse to be unchanged from runs with the same grid resolution and only four particles per 
cell. The overdense seeds of the gravitational instability are under all circumstances resolved 
with several hundreds or thousands of particles, explaining this result. 

Although we have clearly not yet obtained convergence in the gravitational collapse of the 
particles, the consistent trend towards higher peak densities and gravitational instability at 
lower average column density with increasing resolution lends strong support to our hypothesis 
that this mechanism can drive planetesimal formation at gas column densities characteristic 
of the minimum mass solar nebula. 

To compensate for underresolving the peak densities, one could manually sharpen the un- 
derestimated density modes^, but we prefer not to solve for the gravitational acceleration 
better than the pressure term, to avoid any risk of artificial fragmentation^. On the other 
hand, studies that use the adaptive mesh refinement method^ in combination with a kinetic 
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Figure 6: Maximum particle density on the grid, assigned from particle positions with the 
TSC scheme^ 133 , versus time for runs with no self-gravity and no collisional cooling. The 
maximum density goes up by a factor of 3-4 for the 128 3 simulation, in part because the density 
field assigned with the TSC scheme on the grid has less smoothing at higher resolution. 

particle code with collision detections^ could be adapted to fully resolve the gas concentra- 
tion of particles, and then follow the collapse and gravitational fragmentation of the boulder 
layer all the way up to solid density, allowing determination of the minimum column density 
for planetesimal formation and the final multiplicity of the self-gravitating boulder clusters 
that we model. 




1.7 Collisional cooling 

In the runs presented in the main text, we include the effect of collisional cooling on the 
collapse of the mid-plane layer into gravitationally bound clusters. However, the collisional 
cooling only marginally changes the column density at which collapse initially occurs. To 
demonstrate this, we have run simulations without collisional cooling to quantify its impor- 
tance for gravitational collapse (see column 8 of Table [2|. We show in Fig. [7] the maximum 
bulk density of solids in simulations with and without collisional cooling, for 64 3 (left plot) 
and 128 3 (right plot) grid simulations. At 64 3 resolution collapse occurs at G = 0.4 both 
with and without collisional cooling, although its inclusion seems to cause collapse to hap- 
pen somewhat more quickly after self-gravity is turned on at t = 0. At a resolution of 128 3 
zones, collapse occurs in the presence of collisional cooling at G = 0.2, whereas the column 
density limit without collisional cooling is roughly 50% higher. Collisional cooling is thus not 
a prerequisite of the collapse, but rather allows collapse to occur in somewhat lighter discs. 
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Figure 7: Comparison of runs with (bright /yellow) and without (dark/blue) collisional cooling 
(self-gravity and collisional cooling are turned on at t = 0). At 64 3 resolution, we find that 
collapse occurs at the same gas column density, independent of collisional cooling, although 
collisional cooling lets the collapse happen faster. At 128 3 collapse with collisional cooling 
occurs at G = 0.2, while G = 0.3 is needed when collisional cooling is not applied. 



1.8 Jeans and Toomre criterion 

The Jeans criterion^ 51 * 54 ! states that the local Jeans length for gravitational collapse of gas must 
be resolved by at least a few grid points, as otherwise artificial fragmentation may occur in 
marginally stable structures because of numerical discretisation errors. This is particularly a 
problem in the presence of rotation 54 ! 55 ! A corresponding Toomre^ criterion^ is appropriate 
for accretion discs. We show in §1.10 that the initial radial collapse phase is resolved by 



around 20 grid points, but gravitationally bound clusters eventually reach the size of a few 
grid cells, for which the unstable wavelengths that would lead to further fragmentation are 
not resolved. It is not snown whether artificial fragmentation is an issue for solid particles 
moving on a fixed mesh. The TSC scheme inherently smooths out small scale power in 
the assigned particle density^ 3 -^ (see also [1.2), so the contribution of the small scales to 



the gravitational potential is underestimated, whereas the particle velocity dispersion that 
counteracts contraction is resolved equally well at all scales. It is thus more likely that our 
numerical scheme suppresses fragmentation at the small scales and at crude resolutions; this 
is also supported by the resolution tests where gravitational collapse occurs at lower column 
density in higher resolution models. We have no reason to believe that the collapse of the 
particle layer should stop anywhere near the grid scale: the velocity dispersion in the overdense 
clusters is around 1% of the sound speed, so the pressure support is negligible^, and drag force 
and inelastic collisions are highly efficient at dissipating kinetic energy. Instead we consider 
the condensed particle clusters to be equivalent to (numerically expensive) sink particles, a 
good approximation as long as there is no feedback from the unresolved scales back to the 
large scales^. 
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1.9 Gammie cut-off 



Gammie^ argues that one must exclude wave numbers above a limit from the solution to 
the Poisson equation in the shearing box in order to keep the gravitational acceleration ap- 
proximately isotropic at small scales in the presence of shear. The cut-off is at fc max = \/2&Ny 
where k-^ y is the Nyquist wavenumber. This cut-off has been applied in all our simulations. 
We have run test simulations without the cut-off and found no significant differences in the 
resulting gravitational fragmentation, probably because the TSC assignment scheme itself 
underestimates the potential at these small scales. 



1.10 Size of the forming bodies 

The linear stability analysis of Goldreich & Ward- gives the largest wavelength unstable to 
radial self-gravity as 

Agw = o2 P , (44) 

where G is the gravity constant, U p is the column density of solids and i?K is the Keplerian 
orbital frequency. This expression is formally valid in the limit of vanishing particle velocity 
dispersion and ignores the potentially important effect of drag forces on the collapse* 42 * 58 *. 
Inserting nominal values at r = 5 AU in the protosolar nebula gives 

^■^M^HwX^u) 3 ' (45) 



Using a scale height of H = 3 x 10 12 cm (see { 1.5 1 and U p = 15 g cm 2 in the radial overden- 



sities of our simulations yields \qw/H w 0.1, which is initially well-resolved (the grid size is 
5x = 5y = 5z = 0.005-ff in the 256 3 simulation). Assuming that all the solid mass within a 
single unstable wavelength collapses to a gravitationally bound solid object, the radius of the 
object is 

R w ( ^exi w Y /3 w 50km£2/3 (MA-** ( ^ \ (jy\ 2 ( p- y 1/3 

\ p. J \M Q J U-5gcm- 2 y V5AUy V2gcm-V 

(46) 

where £ < 1 is a parametrisation of the most unstable wavelength relative to the largest 
unstable wavelength. A reasonable value is £ = 1/2. 

At r = 1 AU with S p = 15 gem -2 one recovers the 10 km planetesimals found by Goldre- 
ich & Ward . The size of planetesimals falls to less than 1 km if the velocity dispersion 
of the particles is included, but the derivation ignores the effect of drag forces on energy 
dissipation ! 59 * 60 * 4 - 2 ^ 5 ^. Considering, on the other hand, an orbital radius r = 5 AU, with 
U p = 30 g cm -2 in the radial overdensities found in the current work, the planetesimals reach 
radii of around 600 km, because the unstable wavelengths contain far more mass. These sizes 
agree roughly with the results of our numerical models. 
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1.11 Collision speeds 



We show in Fig. [8] the collision speeds between bodies of different sizes in a 128 3 simulation 
without self-gravity or collisional cooling (from a snapshot taken at i = 20T or b). We consider 
a million particles of four different particle sizes at the same time, i?K r f = 0.25, 0.50, 0.75, 1.00 
and two values of the radial pressure support, Av = — 0.05c s and Av = — 0.02c s . 

We have considered collisions between all particles within the same grid cell by measuring the 
relative speed of every single pair of particles, giving (l/2)iV(JV — 1) unique collisions in a grid 
cell with N particles. We note that this may underestimate the actual rate of fast collisions 
since in reality high speed collisions would happen more frequently. We also note that grid 
points with many particles get a high weight in the overall picture, and if these very dense 



grid points have a lower collision speed due to feedback shielding from the gas (see Figs. 13 



and 14), this will also lead to a generally lower collision speed. All collisions were considered 
for simplicity to be head on. The more realistic case of random impact parameters would lead 
to a reduction of the measured collision speeds by more than 50%. We show for comparison 
the collision speeds measured at t = 20T or b for 64 3 and 256 3 simulations in Figs. [9-10 For 
the 64 3 the Mach number of the gas is overall larger, Ma ~ 0.07 compared to Ma ~ 0.05 at 
128 3 , and collision speeds are also larger than in Fig. [8[ There is very little difference between 
the 128 3 and the 256 3 collision speeds, but the data is based on only one snapshot. More 



statistically significant convergence tests, given in £1.11.1, show a slight increase in collision 



speeds (relative to large scale rms speeds) with increasing resolution. 

Equal-sized bodies have no differential drift, so turbulent motions dominate the distribution 
of collision speeds, while differential drift (dotted lines in Fig. [8]) dominates the distribution 
of collisions between bodies of different sizes. 

Following BenzP we calculate the threshold for destructive collisions. We use the specific 
incoming kinetic energy 

K = {l/2)mv 2 /M, (47) 

where v is the collision speed and m and M are the impactor and target mass, respectively, 
with M > m. Collisional fragmentation occurs at a limit JC* ~3x 10 6 ergg -1 for metre-sized 
solid bodies and /C* ~ 6 x 10 5 ergg -1 for metre-sized rubble piles'^. This allows us to calculate 



the destructive collision threshold from equation (47) for any pair of bodies (see Fig. [8[). We 
have ignored any shape effects that could lead to differential drift of equal mass bodies^. 
More collisions exceed the destruction threshold for equal-sized bodies than other cases, but 
destructive collisions still only form a minority of the total, although the fraction of destructive 
collisions depends on poorly known material properties of the boulders. Note that our bodies 
are in fact smaller than the 1 m considered by Benz, so the destruction threshold speeds 
may in reality be higher since smaller bodies are expected to be able to withstand higher 
collision speeds. Higher radial pressure support does not produce markedly more collisional 
fragmentation, as can be seen by comparing the grey and dark lines in Fig. [8] The collisions 
for which the relative speed increases with the drift speed are between different-sized bodies, 
and those collisions can withstand much higher relative speeds, whereas higher radial pressure 
support does not give equal-sized bodies a higher collision speed. A third possibility for the 
internal makeup of the bodies is that they are icy, porous structures. It is shown by Ryan et 
al.^ that such objects are generally as strong as the silicates considered by Benz^. We ignore 
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Figure 8: Distribution of collision speeds between boulders of different sizes for a radial 
pressure support of Aw = — 0.05c s (grey line) and Av = — 0.02c s (dark line), both for 128 3 grid 
points and 1,000,000 particles without either collisional cooling or self-gravity. The laminar 
solution for the relative speed due to differential drift (in radial and azimuthal velocity) is 
shown with a dotted line (with an assumed solids-to-gas ratio of e = 0.5). This is the major 
contribution for different-sized bodies, whereas equal-sized bodies only obtain collision speed 
due to turbulent motions in the gas. The threshold for destruction is shown for solid bodies 
(dot-dot-dot-dashed line) and for rubble piles (dot-dashed line) following BenzP. Equal-sized 
bodies have the lowest destruction threshold, but that is still in the far wing of the collision 
distribution. Small bodies and large bodies can survive collisions at much higher speeds, 
sometimes even beyond the range shown along the 5v axis (i.e. > 0.1c s , where c s =500 m/s). 
The collision speed between equal-sized bodies is relatively unchanged when increasing the 
pressure support to Av = — 0.05c s . The increased differential drift, however, manifests itself 
in higher collision speeds between different-sized bodies. 
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Figure 9: Same as Fig. [8] but for 64 3 grid points with 125,000 particles. The overall Mach 
number of the turbulent flow is 30-40% higher than in the higher resolution runs, hence 
collision speeds are somewhat higher. 
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Figure 10: Same as Fig. [8] but for 256 3 grid points with 8,000,000 particles. Here we consid- 
ered for simplicity only 1/8 of the immense amount of possible collisions. The collision speeds 
agree well with those of Fig. [8} 
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Figure 11: The distribution of particle collision speeds (volume average, left plot) in units of 
sound speed, also shown weighted with number of particles in the grid cell (particle average, 
right plot), for Av = — 0.02c s and four different particle sizes i?K = 0.25, 0.50, 0.75, 1.00. The 
distribution function / is normalised such that f(t, 5v)d(5v/c s ) = 1 for all t. Collisional 
cooling and self-gravity are turned on at t = 0, halfway through these runs, with little effect 
on average collision speeds until late times. Collisional cooling is only efficient in the few 
regions where the solids-to-gas ratio is very high, so the average collision speeds are not much 
affected. Accretion onto the condensed out clusters is visible in the right panel at late times. 



in the present analysis other effects of collisions, such as erosion that may deplete boulders 
of mass even if the collision speed is below the threshold for destruction. 

A quite serious issue related to collisional fragmentation is that it could reduce the number 
of particles large enough to take part in gravitational collapse prior to the collapse to as little 
as 20% of the total mass of the solids'^SES. This is particularly relevant in the inner solar 
nebula where the sound speed is high. In that case we speculate that an augmentation of 
the solids-to-gas ratio by up to a factor 5 would be needed for gravitational collapse to occur 
in the midplane layer. Note that we already assume that half of the solid mass is present in 
small grains too well-coupled to the gas to take place in the collapse (assuming that the total 
abundance of solids and ices is 0.017 times that of the ga^SJ). A more detailed modelling of 
collisional fragmentation is needed to quantify its actual effect on the gravitational collapse. 



We show in Fig. 11 the time evolution of the average collision speed of the particles for a 
radial pressure support of Av = — 0.02c s . The volume average collision speed distribution 
has one speed associated with each grid cell whereas the particle average collision speed 
distribution has the local collision speed weighted by the number of particles in the grid 



cell. Collisional cooling, turned on together with self-gravity at t = in Fig. 11 does not 
reduce collision speeds broadly, since only regions with very high solids-to-gas ratios have 
sufficiently high collision rates for cooling to be effective. Accretion with high relative speeds 
onto condensed clusters is visible at late times in the particle average distribution. Some of 
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the high-speed collisions that occur during the accretion phase, especially among boulders 
that have already been deformed in earlier collision events, may in reality lead to erosion or 
destruction of the colliding bodies 3 . However, the fragments remain within the Hill sphere of 
the gravitationally bound cluster, so their dynamics is dominated by the gravity of the nearby 
cluster rather than by the tidal force of the protostar. Furthermore, the mean free path of 
the fragments, i/H « ]?KTf/(Pp//Og) ~ 10~ 3 , is shorter than the size of the bound clusters, so 
they will undergo collisions with the remaining boulders, leading only to the further growth 
of the remaining boulders as they sweep up the fragments. All in all collisional destruction of 
boulders (or, rather, the lack of widespread collisional destruction) is important for obtaining 
an initial condition with a sufficient amount of the solid material bound in boulders, whereas 
collisional destruction during the gravitational collapse phase is probably only of secondary 
importance. 



1.11.1 Comparison of rms and collision speeds with analytic values 



We now aim to explain why the relative (i.e. collisional) speeds are only 5v ~ 0.01c s for 
marginally coupled particles despite turbulent motions of Mach number A4 = u t /c s = 0.05. 
We summarise the four main effects: (1) due to finite friction times particles are not fully 
accelerated by the turbulence, so rms particle speeds are reduced below that of the gas; (2) 
particle drag on gas tends to reduce both random and relative speeds; (3) particles are well 
enough coupled to eddies that collision speeds are slower than rms speeds due to entrainment 
by the same extended eddies; and (4) even though we have a fairly long inertial range in the 
256 3 simulations, particles gain relative speed from scales near the onset of the dissipative 
subrange, so it is not certain that all the relevant scales are resolved. We present extensive 
documentation of effects (1) and (2) in this section, while point (4) is given some support by 
a 10-20% increase of collision speeds relative to rms speeds when going from 64 3 to 128 3 grid 



points (see Figure 16) 



We emphasise that there is no complete theory describing all these effects, so the current 
numerical study deepens our understanding of these issues. Notably an analytical theory for 
collision speeds that includes orbital dynamics is, to our knowledge, missing (even without 
drag of particles on gas), so verification of point (3) above is difficult and the point is left as 
speculation in this work. We also emphasise that our results show that the collision speed 
increases slightly with increasing resolution and that higher resolution studies will be needed 
to see if the collision speeds have converged. The important issue of collision speeds of 
marginally coupled solids will need more consideration in the future. 

First a word on nomenclature and on calculation procedure: 



• With rms speeds we refer to the overall root-mean-square speed of either gas or particles. 
This value is directly obtained from calculating the root-mean-square of all values of a 
chosen component of gas or particle velocity. 

• Collision speed refers to the local velocity differences of particles within a single grid cell 
(similar to the sound speed or mean thermal speed of the gas). We calculate the mean 
collision speed of particles by following 1000 particles from the greater ensemble and 
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With feedback (64 3 ): 
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Table 3: Particle and gas rms speeds and their ratios, for different values of the friction time. 
The top part of the table shows results for passive particles (i.e. with no friction force on 
the gas), whereas the bottom part shows the result when gas is allowed to feel friction. The 
gas Mach number is around Ma ~ 0.05 for 128 3 grid points, while Ma ~ 0.07 for 64 3 grid 
points. Small particles show the expected trend towards equal particle and gas rms speeds 
(the ratios are shown in the last three columns). The azimuthal particle rms speed falls 
quickly with increasing friction time as the loosely coupled particles go on elliptic orbits, and 
the ratio between radial and azimuthal particle rms speed approaches the expected value of 
0.5. Particle rms speeds fall by a significant factor when gas is allowed to feel the friction 
from the particles - for i?KTf = 1.0 the reduction is around 25% when adding the components 
quadratically. 
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Figure 12: Comparison of simulated particle rms speeds and scale-heights to a simple ana- 
lytic model of turbulent forcing that includes orbital dynamics, but no feedback^. Symbols 
are the measured rms particle speeds from single particle size simulations r s = ^K^f = 
0.2,0.5,1.0,2.0,5.0 (128 3 grid points and 1,000,000 particles) either without (black symbols) 
or with (grey symbols) particle feedback. Left plot: Particle scale-height (crosses) and vertical 
rms speeds (asterisks) compared to analytical curves which use the measured vertical rms 
gas speed, u z = 0.024c s and an eddy time, i e ddy = l.Oi?^ 1 as input. Right plot: Measured 
radial (crosses) and azimuthal (squares) particle rms speeds compared to analytic curves us- 
ing the measured gas velocities (radial, u x = 0.032c s ; azimuthal u y = 0.030c s ; and correlated 
(u x u y ) 1 / 2 = 0.016c s ) and i e dd y = 1-OJ?^ 1 as input. See text for further discussion. 



measuring the velocity components of each particle relative to a random other particle 
in its grid cell, calculating in the end the rms of each velocity component from the 
distribution of relative velocities. We make sure that the particle has a neighbour by 
excluding grid points with only 1 particle. We also exclude a neighbour particle if it is 
one of the 1000 chosen particles. This method is computationally much cheaper than 
the full tabulation of all collisions used in Fig. [8] 



For all measurements presented in Tables [3] and [4] we have averaged over snapshots taken 
equidistantly 10 orbits apart, starting at t = 20T or b where the sedimentary mid-plane layer 
has already had 10 orbits to form in balance with the turbulent diffusion. 

We show in Table [3] the rms speeds of gas and particles for different values of the 
friction time (from 128 3 simulations with single-sized particles). The top part shows the 
results of simulations where drag force from particles on gas was not applied, whereas this 
back-reaction drag force is included in the bottom part of Table [3| The turbulent rms speed 
of the gas has in all cases a Mach number of around Ma = Ut/c s = 0.05 for 128 3 grid points 
(seen by adding up the three components quadratically) , while Ma = 0.07 for 64 3 grid points. 
The smallest particles (fi-aji = 0.2) have comparable rms speeds, but the rms speeds fall with 



stopping time as expected (see Fig. 12 and explanation below). 
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Figure 13: Gas rms speed as a function of height over the mid-plane in simulations with 128 3 
grid points and 1,000,000 particles. At each z- value we have calculated the dispersion of the 
three velocity components. There is a clear dip in the gas velocity dispersion around the 
mid-plane where the solids-to-gas ratio is high. This dip is not present in the run with no 
back-reaction friction force (lower right panel). 

To understand the simulated rms speeds of particles in MRI turbulence, we compare to Youdin 
& Lithwick (hereafter YL)^. YL consider the rms response of particles to stochastic turbulent 
forcing (with a temporal Kolmogorov spectrum) with orbital dynamics - in particular epicyclic 
motions - included, which previous works have ignored . The inputs for the simple analytic 
theory include the measured rms speeds of the turbulence (radial, azimuthal, vertical and 
the radial-azimuthal correlations), the friction time of the particles, and the eddy time, i.e. 
the correlation time of turbulent velocity fluctuations. Since the eddy time is a measured 
quantity, we treat it as a free parameter and find the reasonable value t e ddy ~ • 

The left panel of Fig. [12] shows that vertical speeds and scale heights of particles agree very 
well with the analytic theory of YL. Since vertical particle motions are decoupled from in- 
plane motions, the analytic description is fairly simple. When feedback is included (grey 
symbols) the random particle velocities drop, by around 25% for r s = 1.0 particles. This is 
not due to a decrease in the total turbulent Mach number (which actually goes up slightly) , 
but likely comes from the increased particle inertia in dense clumps. 
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Figure 14: Normalised Fourier amplitude of the gas velocity as a function of wavenumber k 
(in units if the largest scale in the box ho) and height over the mid-plane z. Shown here is 
a simulation with single-sized J?kTf = 1 particles (no self-gravity or collisional cooling). We 
have done 1-D Fourier transforms along the x-direction and averaged this over all values of 
azimuthal coordinate y and time t for each given height z. The power spectrum has been 
normalised with the average value at each wavenumber k. The colours represent values of 
±20%. Feedback drag force from particles on the gas reduces the velocity amplitude in the 
mid-plane by around 20%, explaining a similar reduction in rms speeds visible in Table [3} 



The right panel of Fig. 12 shows in-plane particle velocities compared to the YL model. The 
agreement is good, if not as perfect as the vertical case, due to more complex dynamics. 
Notably, the prediction that azimuthal velocities fall before radial speeds (as r s increases) 
due to epicyclic motion is confirmed. Drag feedback (grey symbols) again produces a decrease 
in rms speeds, by about 30% in the radial speeds for r s = 1.0. 

It may be surprising that the eddy time-scale is so short, t e ddy ~ I^k* ( a sim i iar result for 
the eddy time was found recently by Johansen, Klahr, & Mee^ by considering the relation 
between turbulent diffusion and gas rms speeds). Simple mixing length estimates based on 
Fig. [3] yield an eddy time at least an order of magnitude higher than our best fit value, 
i.e. tfc = l/(kuk) ~ lOO,!?^ 1 on the largest scales. But the large scales of the flow are 
dominated by Coriolis force and centrifugal terms rather than non-linear advection, and thus 
the correlation time of these important structures is forced to be similar to the Keplerian 
frequency* 64 * 18 l Mixing length theory gives a good estimate of the life times of large scale 
structures, but their correlation time, which is the important quantity for mixing purposes, 
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Figure 15: The root-mean-square of the particle z-coordinate as a function of time, with and 
without feedback, for a 128 3 simulation with single-sized J?KTf = 1 particles. The scale height 



of the particles falls by around 20% when feedback is applied, in accordance with Fig. 14 



that shows that feedback damps MRI turbulence by around 20% in the sedimented mid-plane 



layer (the reduction in scale height is also seen in the left panel of Fig. 12). 



is never longer than a Keplerian shear time i? K . 

The net effect of decoupling and feedback gives a total (radial, azimuthal and vertical) rms 
speed of 0.027c s for r s = 1.0, nearly half the gas rms speed of vt = 0.05c s . Regions of high 
particle density damp out the gas velocity on the time-scale r g = Tf/(1 + e) where Tf is the 
friction time of the particles and e is the solids-to-gas ratio. It was shown by Dobrovolskis 
et al.^H that this damping can lead to a minor reduction of turbulence around the mid-plane 
layer. We show in Fig.[l3]the gas velocity dispersion as a function of height over the mid-plane. 
There is a clear, though rather modest, dip around the mid-plane where the solids-to-gas ratio 
is also high. The dip is only visible in simulations that include the back-reaction friction force. 



In Fig. 14 we show the magnitude of the different Fourier components of the gas velocity as a 
function of height over the mid-plane. Here we have Fourier-transformed along the x-direction 
and averaged over every single y coordinate at each height z. The spectrum has afterwards 
been normalised with the average spectrum so that differences appear more clearly. There 
is again a clear dip, of around 20%, in the strength of the velocity fluctuations around the 
mid-plane. The dip occurs at all scales, although the small scales of u v are less affected than 



the large scales. In Fig. 15 we show the particle rms z-coordinate as a function of time. With 
particle feedback on the gas there is a drop by around 20% in the scale height of the particles, 
in agreement with what is expected for a 1 — 0.2 2 = 36% drop in the diffusion coefficient^. 
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Without feedback (64 3 ): 
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With feedback (128 3 ): 
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Table 4: Particle rms speeds (a^, ), collision speeds {5v x \ y \ z ) and their ratios. The top part of 
the table shows results without drag force on the gas, the bottom part includes drag force on 
the gas. The collision speeds are very low for small particles, but the ratio between collision 
speed and rms speed increases with increasing friction time as the equal-sized marginally 
coupled particles are big enough to have significantly different histories when they meet. 
Back-reaction friction force on the gas has very little influence on the ratio of collision speeds 
to rms speeds. There is a 10-20% increase in the ratio of collision to rms speed when going 
from 64 3 to 128 3 grid points, possibly due an expansion of the inertial range of the turbulence. 
Higher resolution studies will be needed to determine whether this trend continues of whether 
the collision speeds have already converged. 
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Figure 16: The ratio of particle collision speed 5v p to large scale particle rms speed a p , both 
calculated as quadratic sums of the directional components given in Table [4] and shown for 
two resolutions and with and without particle feedback on the gas. The collision speeds 
increase as expected with increasing friction time as particles decouple from fast large scale 
eddies. Collision speeds increase by 10-20% when increasing resolution from 64 3 to 128 3 , 
likely because the inertial range of the gas turbulence expands to smaller scales. Feedback 
has only little influence on the ratio of collision speed to rms speed, but the overall rms speed 



of particles is reduced by feedback (see Fig. 12). Higher resolution studies will be needed to 
determine whether J?KTf = 1 particles will eventually obtain the fully mixed dv p /a p = \[2. 



1.11.2 Why are collision speeds so low? 



Table [4] shows the collision speed 5v of the particles for comparison with the particle rms 
speeds. Here we have chosen a thousand particles throughout the simulation and calculated 
their relative speed with a random different particle in the same grid cell. The collision speed 
of marginally coupled particles (with i?KTf = 1) is around 50% of the overall rms speed. This 
factor 2 decrease in the relative speeds on small scales compared to rms speeds remains to be 
explained. An exact analytic prediction of this effect is difficult, and has not to our knowledge 
been done with either orbital dynamics or feedback, let alone both effects. Higher resolution 
simulations will be needed to determine if small scale eddies can change this result. Fig. [3] 
indicates that over an order of magnitude in the inertial range is resolved, but scales around 
the onset of the dissipative subrange may be important for collision speeds. 

We can think of two possible reasons for the low collision speeds. The first is insufficient 
numerical resolution. The numerical dissipation of the code results in an underestimate of 
the amplitude of turbulent velocities at scales smaller than around 6 grid points. At the same 
scales, the particle-mesh drag force (£1.2) begins to be underestimated^. Taking the stopping 
length of the boulders to be £ st0 p = Tft> rms , where v Tms potentially depends on Tf, this yields for 
marginally coupled particles a stopping length of £ s t Q p = 0.025-ff, when reduction of the rms 
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Figure 17: Particles in the densest grid cell at t = 20T or b for the 256 3 simulation presented in 
the main paper (this is before collisional cooling and self-gravity are turned on). The panels 
show a gradual zoom out on the surrounding grid cells (still in the same z-plane) . No spurious 
structure is evident either at the subgrid scale or at the interfaces between grid cells, since the 
high order particle-mesh drag force scheme samples 27 nearby grid points when evaluating 
the gas velocity at the position of a particle. The particles are shown in velocity space in the 



left panel of Fig. 19 



speed due to marginal coupling and feedback is taken into account. The resolution is either 
5x = 0.01, for 128 3 grid cells, or Sx = 0.005 for 256 3 grid cells. If we assume that relative 
collision speed between particles of stopping length £ stop is typically given by turbulent eddies 
of wavelength 2£ s t op , then these eddies are of size 5-10 grid cells, depending on resolution. 
Thus there is potentially a damping of the collision speeds, as some of the scales that are 
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Figure 18: The mean particle density in the mid-plane (left panel) for a simulation with 
single-sized J?K T f = 1 particles and feedback on the gas. The right panel shows the radial 
velocity w x of the particles in the mid-plane. Isolated particles would have w x = — 0.02c s , but 
the increased inertia in the mid-plane has decreased the radial drift by a factor two there. The 
dashed line shows the expected radial velocity calculated from the mid-plane density shown in 
the left plot. There is relatively good agreement, an indication that drag forces are calculated 
correctly in the code (we note that a thorough test of the particle-mesh drag force scheme was 
performed by Youdin & Johansen® ). Radial drift is a bit lower than expected, but that may 
be explained by large scale gas density fluctuations that modify the radial pressure gradient 
enough to slow the radial drift of marginally coupled particles^. 



important for relative motion are just at the onset of the dissipative subrange. Table [4] and 



Fig. 16 show that the ratio of collision speeds to rms speeds does increase by 10-20% when 
going from 64 3 to 128 3 grid points, although the increase is much more modest when feedback 
is included. An increase in grid resolution, implying both an increased inertial range and a 
sharpening of the particle-mesh scheme at small scales, will be needed to determine whether 
the collision speeds will continue to go up with increasing resolution. 

A second possible explanation is missing orbital dynamics in analytical estimates. The colli- 
sion speeds predicted by Voelk et alJ™ rely on a Kolmogorov spectrum with eddy turn over 
times from mixing length theory. This approach does not take into account epicyclic motion, 
nor the upper limit to structure correlation times given by the Coriolis force. YL already 
showed that orbital dynamics are important for particle rms speeds, but there is no similar 
theory for collision speeds that takes into account orbital dynamics. Such analytical work 
should be a high priority for future research because of its importance to the dynamics and 
growth of boulders across the m-sized barrier. 

We note that the particle-mesh drag force scheme of the Pencil Code has been extensively 
tested in the recent paper by Youdin & Johansen^. The same TSC assignment /interpolation 



scheme is also used for self-gravity, for which we present a test problem in { 1.3.1 
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Figure 19: Particles shown in velocity space (v x ,v y ). The left panel shows the particles from 



Fig. 17 - the colours refer to the 16 individual grid points in the lower right panel of Fig. 
17 The right panel shows a 128 3 simulation with single-sized /?K T f = 1 particles. Notice 
the higher axis range. Here we have chosen 20 random grid points containing between 500 
and 1000 particles. The overall spread of particle velocities is large, comparable to 0.02c s in 
each direction, whereas the spread within individual grid cells is much smaller. The positive 
correlation between the two velocity components arises from the positive Reynolds stress of 
the MRI turbulence. 



We proceed with a few reality checks of particle positions in physical space and in velocity 
space. While we do not claim that any of these tests are exhaustive, they are meant to serve 
the purpose to exclude that spurious particle structure at the grid scale is the source of the 
low collision speeds. We show in Fig. [17] particle positions around the densest grid point at 
t = 20T or b for the 256 3 simulation presented in the main paper, before collisional cooling 
and self-gravity are turned on. No spurious structure is evident at subgrid scales or at the 



interfaces between the grid cells. In Fig. 18 we show the mean particle density in the mid- 



plane as a function of time (left panel), for a simulation with single-sized J?K T f = 1 particles. 



The right panel of Fig. 18 shows the mean radial drift velocity in the mid-plane. Particles 
drift slower than the expected — 0.02c s due to increased particle inertia in the mid-plane^. 
The dashed line indicates the expected drift velocity given the mid-plane density shown in 
the left panel. It is in relatively good agreement with the measured drift, another indication 
that the code applies drag forces correctly. A slight decrease in the drift speed may be due 
to large scale pressure gradients in the gas slowing down the radial driftP. 

Finally we show in Fig. [19] particle positions in velocity space (v x , v y ). While the overall spread 
in velocities is large, the individual grid cells (shown with individual colours) show a much 
more coherent velocity structure. We hope to explain this difference in future publications, 
both with increased numerical resolution and with improved analytical models that include 
orbital dynamics on the calculation of the collision speeds. 



36 



1.12 Low ionisation discs 



A high enough ionisation fraction is needed for the magnetorotational instability to oper- 
ated This threshold may be reached in hot parts of the disc and in parts where the column 
density is low enough that cosmic rays penetrat d 68 * 69 * 70 *. Typically a penetration depth of 
Eg ~ 100 g cm" 2 is found. Other sources of ionisation are decay of radioactive elements and 
chemical reactions. Small dust grains, on the other hand, soak up electrons and reduce the 
conductivity^. This has lead to the concept of a dead zone in the midplane of the disc where 
ionisation is insufficient for MRI 2 -D. Several mechanisms have been proposed for reviving the 
dead zone: growth of the peak dust grain size by an order of magnitude^, using the turbulent 
kinetic energy to ionise molecules, thus maintaining any turbulence that is already present^, 
and X-ray flares from the central star that would ionise the disc periodically^. 

We investigate the limit in which the MRI completely fails to influence the midplane by 
modelling the same problem as in the main text, but absent magnetic fields. A more advanced 
model would still have magnetically active surface layers that could influence the "dead zone" 
dynamically* 74 * 75 -^. We do not, in general, expect density fluctuations in the gas similar to 
what happens in the MRI. This is because the Mach number of Kelvin- Helmholtz turbulence 
and streaming turbulence is generally smaller than for the MRjti!ll2l This situation may 
change in a more advanced model with active surface layers that can influence the magnetically 
dead mid-plane, as the surface layers send density waves through the dead zon d 74 * 7 § *. 



1.12.1 Low radial pressure support 



We show in Fig. 20 the evolution of the sedimented mid-plane layer for the usual range 
of particle sizes and a radial pressure support of Av = — 0.02c s . The box size was set to 
L x = L y = L z = 0.1H, 1.32 times smaller than in the magnetic runs, in order to resolve 
the thin mid-plane layer that forms in the absence of global MRI turbulence. A turbulent 
state develops anyway since vertical shear in the radial and azimuthal gas flow is unstable 
to the Kelvin-Helmholtz instability* 5 * 77 * 7 * 13 * 78 * 11 ^. The streaming instability is also a source of 
non-linear dynamics and turbulenc d 10 * 37 * 1 ^*. 



It is evident from Fig. 20 that the particles quickly gather at a single radial location and drift 
radially together for the full duration of the simulation (50 orbits, the Courant time-step in 
these small boxes limits the possible simulation times compared to the magnetic runs). This 
can be understood as an effect of the streaming instability. Any isolated particle drifts quickly 
into the overdense band where the dynamics is so dominated by particles that the radial drift 
is significantly reduced. Very high particle densities are reached within the band, with peaks 
at more than three orders of magnitude times the gas density. 

The high particle densities are very susceptible to gravitational collapse. We show in Table [5] 
that gravitational collapse is possible in discs with as low mass as one half of the MMSN. 
The accretion rate of the most massive bound cluster is generally 3 orders of magnitude lower 



than in the magnetic runs, which follows the expected trend of equation ( 46 ) that the mass 
of the bound object scales with the cube of the column density. The solid size of the most 
massive bound cluster in the 128 3 run of Table [5] is 150 km, smaller than in the runs with 
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Figure 20: Particle column density for 64 3 and 128 3 simulation with no magnetic fields and 
the usual range of boulder sizes. Particles quickly end up in a single azimuthally extended 
roll and slowly drift in together for at least 50 orbits. The maximum bulk density of particles 
reaches 1000 times the ambient gas density. 



magnetorotational turbulence presented in the main text, but still larger than the classical 
view of planetesimals. 



1.12.2 Moderate radial pressure support 



We show in Fig. 21 a comparison between two 3-D simulations: one with magnetic fields 
and one without. Both have a resolution of 128 3 grid points, marginally coupled particles 
with J?K r f = 1 an d a moderate radial pressure support of Av = — 0.05c s . The box size is 
L x = L y = L z = 1.32H for the MRI run and L x = L y = L z = 0.2H for the run with no 
magnetic fields. After an initial period where strong density enhancements form, the solids- 
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Table 5: Resolution study for runs with no magnetic fields and Av = — 0.02c s . Col. (1): Mesh 
resolution. Col. (2): Number of superparticles in millions. Col. (3): Number of orbits with 
self-gravity. Col. (4): Minimum self-gravity parameter where gravitationally bound clusters 
form (MMSN has G « 0.05 at r = 5 AU). Col. (5): Corresponding Toomre Q m 1.6G -1 . Col. 
(6): Number of clusters at the end of the simulation. Col. (7): Accretion rate of the most 
massive cluster in Ceres masses per orbit. 



MRI(128 3 ) NoMRI(128 3 ) 128 3 




Figure 21: Comparison of azimuthally averaged boulder column density for models with 
magnetic fields (left plot) and without magnetic fields (middle plot) for 128 3 grid points and 
marginally coupled particles with i?K r f = 1- The model with no magnetic fields initially 
shows some clumping and overdensities, but these die out after 20 orbits after which the 
maximum solids-to-gas ratio stays below 10 for the duration of the simulation with no signs 
of the high density events that are obvious in the MRI run. 



to-gas ratio never rises above 10 in the simulation without magnetic fields. In Fig. 22 we show 
the time evolution of particle density averaged over the azimuthal direction for the moderate 



pressure support non-magnetic run presented in the middle panel of Fig. 21 The particle 



density reaches a configuration where a single standing wave dominates the box, with no 
significant overdensities. This state may come about because of the streaming instability. 
It was observed by Johansen & Youdin^ that marginal coupling produces a very strongly 
turbulent state (with vertical diffusion coefficients similar to those of the MRI). This state 
was nevertheless very clumpy in the simulations of Johansen & Youdin^, but the inclusion of 
vertical gravity in the current work may mean that the correlation times are so short (due to 
vertical oscillations) that clumping is suppressed. The simulations by Johansen & Youdin^ 
also showed that a higher background solids-to-gas ratio leads to a dramatic increase in growth 
rate of the streaming instability, and a decrease in the unstable wavelengths. We show in the 
next section (j |1.12.3 1 that a slight increase in the global solids-to-gas ratio indeed leads to 
significant particle clumping. 
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Table 6: Resolution study for runs with no magnetic fields, Av = — 0.05c s and eo = 0.03. 
Col. (1): Mesh resolution. Col. (2): Number of superparticles in millions. Col. (3): Number 
of orbits with self-gravity. Col. (4): Minimum self-gravity parameter where gravitationally 
bound clusters form (MMSN has G ~ 0.05 at r = 5 AU). Col. (5): Corresponding Toomre 
Q ~ 1.6G . Col. (6): Number of clusters at the end of the simulation. Col. (7): Accretion 
rate of the most massive cluster in Ceres masses per orbit. 



We have tested the effect of self-gravity in the 3-D non-magnetic models with moderate radial 
pressure support as well, with the four different particle sizes we used in the models described 
in the main letter. We have not been able to form accreting self-gravitating clusters, like 
those we observe in discs with MRI turbulence, for any reasonable column density G < 1 at 
either 64 3 or 128 3 resolution. Global magnetorotational turbulence thus appears to overall 
promote rather than impede gravitational collapse of marginally coupled boulders. 



1.12.3 Higher solids-to-gas ratio 



We performed simulations with increasingly higher solids-to-gas ratio with the expectancy 
that the turbulent state would get weaker and overdensities larger. This behaviour was 
indeed confirmed. In Fig. 23 we show the evolution of the mid-plane layer for moderate radial 
pressure support and a background solids-to-gas ratio of 0.03. Extremely high overdensities, 
up to three orders of magnitude higher than the gas, now occur. We show in Table [6] the 
gravitational collapse for the moderate pressure support case. The accretion rate is around 
six times higher than in the case of low radial pressure support (see Table [5]) , but that is 
probably due to the increased global solids-to-gas ratio in the moderate pressure support 
case. 



One possible way to enhance the global solids-to-gas ratio is by radial drift augmenta- 
tionHSEES This would be especially likely in the case of a magnetically active and warm 
outer disc with a strong radial drift. Drifting into the dead inner regions of the disc the 
boulder velocity field would converge at the rim of the dead zone, leading to particle pileups 
there. Another effect, which to our knowledge has not yet been explored, is that the accret- 
ing surface layers must lose all their solids to the dead zone. Considering an active layer 
with a width of one gas scale height H and an orbital distance r, the probability for a small 
dust grains to accrete all the way to the star without diffusing into the dead zone is only 
(H/r) 2 <C 1. Once a grain is lost into the dead zone it is very unlikely that it will diffuse back 
into the active layers because the turbulent diffusion in the dead zone is assumed to be very 
low. A third possibility for augmenting the solids-to-gas ratio is by photoevaporation of the 
gaseous part of the disc^ 
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Figure 22: Time evolution of the boulder density averaged over the azimuthal y-direction in 
a 3-D simulation with 128 3 grid points and 10 6 marginally coupled particles with /?K T f = 1- 
A Kelvin-Helmholtz instability develops from the initial sedimentation, driven by the shear 
in the radial drift velocity of the gas^H The whole particle layer eventually participates in a 
single standing wave where particles drift and oscillate around the midplane (see main text 
for details). 

1.13 Heating and cooling 

Drag and collisions dissipate the kinetic energy of the solids. In the extreme case, all this 
heat is released locally in the gas. We consider here the case of ineffective radiative cooling 
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Figure 23: Particle column density for a global solids-to-gas ration of 0.03. The particle 
density peaks at more than three orders of magnitude higher than the gas density. 



and an ideal gas equation of state, rather than the isothermal equation of state assumed in 
the rest of the paper. The only means to transport energy in this model is by turbulent heat 
conduction. In reality collisional cooling may transfer most of the energy into deformation of 
the colliding bodies, but we consider here the effect of the most extreme case of local heating 
on gravitational collapse. 

We have ignored viscous heating under the assumption that the gas background has found an 
equilibrium between radiative loses and viscous heating. The extra gas shear and compression 
in the collapsing particle clusters is negligible due to the pressure support of the gas, so any 
extra viscous heating in regions of collapse is not important. We have found no traces of 
the gravitationally contracting clusters in the gas density or in the gas velocity field (gas is 
allowed to exert and feel self-gravity in the simulations that include heating). 



We show in Fig. 24 the gas temperature as a function of time for a run with 64 grid points 
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Figure 24: The maximum and mean temperature of the gas as a function of time in our 
adiabatic model, for models with 64 3 zones and G = 0.5 (left), and 128 3 zones and G = 0.2. 
Self-gravity is turned on at t = 0; the first collapse, seen in the fall in the mean potential of 
the particles, begins after approximately 10 orbits for the 64 3 run. Peaks in the maximum 
temperature coincide well with times where gravitational energy is released (the peak at 
t = —20 is due to compression in the initial random velocity field). The energy is transferred 
to the gas by drag and inelastic collisions. The mean temperature rises negligibly because 
gravitational collapse occurs in only a small fraction of the total gas volume. 



and G = 0.5 (a 25% higher column density was needed for collapse in this case, which we 
attribute simply to the stochastic nature of the collapse) and for a run with 128 3 and the 
usual G = 0.2. Self-gravity is turned on at a time t = 0. We let gas exert and feel self-gravity 
in the simulations with heating, although this causes no real difference in the behaviour. This 



effect is otherwise not included in the simulations (see £1.5.3). The maximum temperature 
rises to no more than 30% over the initial value at times of gravitational collapse (seen in the 
mean gravitational potential energy of the particles). Peaks in the maximum temperature 
coincide well with gravitational energy release events. The peak at t = —20 happens due to 
compression in the initially random velocity field. 



We show in Fig. 25 the dissipation rate of drag and collisions. Drag dominates the volume- 
averaged energy budget, with a heating rate at least six orders of magnitude higher, because 
collisions are only important in the few grid points where the solids-to-gas ratio greatly exceeds 
unity. Assuming a temperature of around 80 K in the non-heated state, the temperature 
reaches a maximum of approximately 104 K in the 64 3 simulation, still too low to even melt 
the icy component of the boulders. 



1.14 Varying radial pressure support 

Cuzzi et al.^ gives a detailed review of the values expected for the radial pressure sup- 
port parameter. We have assumed throughout this work a radial pressure support of either 
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G=0.5 (64 3 ) 



G=0.2 (128 3 ) 








Figure 25: Volume-averaged heating rate by drag (dark/blue) and by inelastic collisions 



(bright /khaki) in the same runs shown in Fig. 24 Drag force dissipates six orders of magnitude 
more kinetic energy than collisions do, although collisional dissipation dominates in regions 
where the bulk density of solids is at least two orders of magnitude higher than the gas (see 
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Table 7: Resolution study for Av = — 0.05c s . Col. (1): Mesh resolution. Col. (2): Number of 
superparticles in millions. Col. (3): Number of orbits with self-gravity. Col. (4): Measured 
turbulent viscosity. Col. (5): Minimum self-gravity parameter where gravitationally bound 
clusters form (MMSN has G ~ 0.05 at r = 5 AU). Col. (6): Corresponding Toomre Q ~ 
1.6G . Col. (7): Number of clusters at the end of the simulation. Col. (8): Accretion rate 
of the most massive cluster in Ceres masses per orbit. 



(dlnP/dlnr)(H/r) = -0.04 or (d In P/d In r) {H/r) = -0.1 (see §L5A\ . The value of H/r 
directly yields the local temperature when folded with the mass of the central object and 
the orbital radius in the disc. We have run test simulations of colder discs as well, with 
H/r = 0.02, corresponding to a four times lower temperature at a given location than models 
with H/r = 0.04. In Fig. 26 we compare the evolution of the maximum particle density 
between the standard run with H/r = 0.04 and the cold run with H/r = 0.02. The two 
curves are statistically indistinguishable. Lowering the disc aspect ratio reduces radial drift 
and decreases the unstable wavelengths of the streaming instability, but this apparently does 
not change the evolution of the solids, at least for moderate changes in temperature. 
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Figure 26: Maximum particle density versus time for standard run with H/r = 0.04 (blue) 
and a colder disc with H/r = 0.02 (yellow). The curves appear statistically indistinguishable. 
A lower disc aspect ratio decreases the radial drift of the solids, while making the unstable 
wave lengths of the streaming instability smaller, but apparently neither of these effects 
significantly changes the particle overdensities in the considered cases. 
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Figure 27: Comparison of the topography of the particle layer for the standard run with 
Av = — 0.02c s (left) with a run with 2.5 times stronger radial pressure support Av = — 0.05c s 
(middle), both at 64 3 grid points and 125,000 particles. The increased radial drift is clear 
in the tilted bands, but the streaming instability still produces strong overdensities in the 
particle layer. The right plot shows the maximum particle density as a function of time - the 
peaks in the particle density are reduced by around 30% in the case of the stronger radial 
pressure support. 

A higher radial pressure support can occur for a given value of H/r if the radial pressure 
support \d1nP/d]n.r\ exceeds unity. We next examine the consequences of a higher radial 
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Table 8: Estimated time-scales, in units of orbits, for collisional fragmentation and for coag- 
ulation in three different states of the mid-plane layer. 



pressure support on the maximum density of solids. We show in Fig. 27 the topography of 
the sedimented particle layer and the maximum particle density of two 64 3 simulations: the 
standard run with (<9 In P/d In r)H/r = —0.04 and a run with higher radial pressure support 
and (d In P/d In r)H/r = —0.1. The same comparison at 128 3 is shown in the main text. 
The increased radial pressure support, Av = — 0.05c s , is very clear in the latter case from the 
tilted particle bands, but the streaming instability still produces radial over densities. We have 
checked the column density required for gravitational collapse in models with At; = — 0.05c s 
and found that approximately a factor two higher value of the column density of gas is needed 
for the collapse (see Table [7J compared to the standard model with Av = — 0.02c s . 



1.15 Other processes 

In order to isolate the effect of self-gravity on the dynamics of the solids we have ignored two 
potentially important collisional effects: coagulation and collisional fragmentation. In this 
section we estimate time-scales and effects of these physical phenomena. 

The time-scale of collisions is 

Tf 

r coii = 7 — ^7 — , — 7 , (48) 

(Cp/Cs)(/9p/Pg) 

where c p is the collision speed of the boulders and Tf is the friction time. Assuming spherical 
boulders gives the relation 

± = 3™- (49) 
a, m, 

between the time-scale for radius doubling and the time-scale for mass doubling, which in turn 
is just the collisional time-scale. We consider regions with solid-to-gas ratios characteristic 
of the average midplane in a turbulent disc, the peak density in the turbulent midplane, and 
an intermediate density in a gravitationally contracting cluster, as described in Table [8} The 
same table gives the time-scales for collisional fragmentation and for coagulation for these 
different regions. 

The average state has a coagulation time of close to 100 orbits, comparable to the time-scale 
for radial drift in. To cross the drift barrier effectively a few coagulation times are needed, so 
the growing solids are lost into the inner disc before they can decouple from the gas. A more 
advanced model would take into account the sweep up of small grains as the boulders drift, 
but essentially the same result is reached^. 

The size of the boulders could increase on the time-scale of only one orbit in the overdense 
particle clusters that occur in transient high pressure regions and due to the streaming insta- 
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bility. Coagulation could go even faster in the gravitationally contracting regions. This may 
be how planetesimals eventually form: by coagulation in an overdense environment set by the 
self-gravity of the solids. The strength of the scenario presented in the main text, though, is 
that it does not rely on coagulation being efficient in same way as pure coagulation models 
dd 28 * 27 l Coagulation is needed to build up marginally coupled solids — metre-sized boulders at 
5 AU or centimetre-sized pebbles at 40 AU — but does not directly drive gravitational collapse. 

Collisional fragmentation may also be prevalent in gravitationally contracting clusters. In 
that case one can speculate that the mode of radius growth is by sweeping up of collisional 
fragments by the few lucky boulders that avoid catastrophic collisions with equal-sized bodies. 
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